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ABSTRACT 

A method for solving the viscous hypersonic flow field around realistic config- 
urations is presented. The numerical procedure for generating the required finite 
difference grid and the two-factored implicit flow solver are described. Results are 
presented for the shuttle orbiter and a generic wing-body configuration at hyper- 
sonic Mach numbers. 


1. INTRODUCTION 

With the continuous development of more powerful computers, computational 
fluid dynamics (CFD) has emerged as a viable tool in understanding complicated 
three-dimensional fluid dynamics phenomena in subsonic, transonic and supersonic 
flow regimes. Hence CFD is now being routinely applied to practical problems to 
supplement wind tunnel experimentation in the design process. The recent renewed 
international interest in hypersonic transport is partly attributed to the confidence 
in using CFD solely in design because of the impossibility of simulating the wide 
range of the actual hypersonic flight conditions in a wind tunnel. This of course 
necessitates the development of robust hypersonic CFD codes which are reliable for 
providing important quantities such as the aerodynamic stresses, heat transfer and 
inlet flow conditions to the engine. 

In this paper we describe and apply a numerical technique which is capable of 
solving the hypersonic flow around three dimensional lifting configurations. In this 
technique, the computational grid is created using an efficient numerical hyper- 
bolic grid generation procedure ( Ref. [ 1 ] ) and the equations of motion (represented 
by the Reynolds averaged Navier-Stokes equations) are integrated using the time- 
dependent factored procedure described in Ref. [2]. 

The present numerical procedure is used to compute the viscous flow around 
the shuttle orbiter at Mach number of 7.9 and angle of attack of 25 degrees and 
around a generic wing-body configuration at Mach number 25 and 5 degrees angle 



of attack. In both cases the fluid was assumed to behave as a perfect gas. The 
calculations were performed on the CRAY-2 to allow the use of fine grids in order 
to resolve the details of the complicated flow field associated with the geometries 
under consideration. 

In the following sections, we briefly describe the grid generation procedure, the 
flow simulation procedure, and the computed results. 


2. GRID GENERATION PROCEDURE 

In computational fluid dynamics, a three dimensional body fitted mesh is highly 
desired and can be generated numerically using the elliptic (Refs. [3, 4]), hyper- 
bolic(Refs.[l,5j), parabolic or hybrid (Refs. [6, 7]) schemes. All these schemes involve 
solving a set of partial differential equations and they all possess the desirable fea- 
ture of allowing the user to specify an arbitrary point distribution along the body 
surface. 

In the present study, the hyperbolic scheme was chosen because it automatically 
creates a body normal grid which is important for viscous flow computations. Also, 
the hyperbolic grid generation method is efficient in terms of computer time and 
memory and the resulting grid can be made orthogonal or close to orthogonal every 
where by proper surface grid specification. The major drawback of the hyperbolic 
scheme is its inability to prescribe an exact location for the outer boundary. How- 
ever, with a minor modification, the user can specify an approximate location for 
the outer boundary as will be discussed later. This limited control on the outer 
boundary shape suffices for a wide range of external flow problems. 

The Hyperbolic grid is generated by solving for the cartesian coordinates x,y,and 
2, in terms of the generalized coordinates and <; where, in general: 

£ = C(x,y,z) 

V = r]{x,y,z) 

here £ is the streamwise coordinate, r) is the circumferential coordinate and f is the 
outward coordinate. 

A set of three equations are required to solve for the three unknowns x,y and 
2. Following Ref.[l], the first two equations are obtained by enforcing orthogonal- 
ity relations between f and f lines and between 77 and f lines, whereas the third 
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equation is obtained from the cell volume specification. The following set of partial 
differential equations result: 


r V^= o ( la ) 

= 0 {lb) 

>»?.?) 

where r is defined as (x,y, z) T 

The above nonlinear system of equations are locally linearized about the previous 
step which results in 


- J~ l = AV 


(lc) 



- n) + Bi6r,{fi + 1 - f;) + = yj+i (2) 


where Ai,Bi and C; are 3x3 coefficient matrices that are formed from the local 
linearization of Eq.(l). The vector gi +i contains the user specified cell volumes 


9l + 1 


0 
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In Eq.(2), 6 is the central difference operator used in f and rj directions whereas V 
is the backward differencing used in f direction. Using approximate factorization, 
Eq.(2), is reduced to 

£r\L${n+\ - n) = c, _1 <fj+i (3) 


where and £c are block tridiagonal matrices given by: 


£„ ={I + Ci~ 1 B l S n ) 


Cc-^il + Cr'AtSc) 

Therefore, the vector r } + 1 is obtained by solving the following sequence: 

£t] 9 l+i = C/ 0/+i 
C^^ri+i = fi + 1 
fi + 1 = ft + Vffi+i 

Starting from the body surface, the mesh is obtained by marching outward using 
Eq.(3). In order to improve stability, numerical dissipation terms are added in the 
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£ and 77 directions. The dissipation consists of second-order implicit and fourth- 
order explicit terms which are scaled to the local mesh spacing ( Ref. [ 1 ] ) - Additional 
optional implicit smoothing(Ref.[8]) is used for complicated body shapes to prevent 
grid line intersection near highly concave sections. 

Depending on the boundary conditions specified at the £ direction, the user can 
create ’’C-O” (Fig. la) or ”H- 0 ” (Fig. lb) meshes suitable for supersonic external 
flow problems. The boundary condition procedure developed in Ref.jljis employed 
at £ — 0 and £ = £ maI • This procedure takes advantage of the fact that the 
eigenvalues of either Cf l A or C^ l B are of the form (cr, 0 ,-ct) where a is real. The 
zero eigen value permits the specification of a combination of the dependent vari- 
ables^. g. z = f(x,y)) or the use of one sided differencing at the boundary. Hence, 
either one of the following two stable boundary procedure can be used. The first is 
to specify one combination of the dependent variables and use one sided differencing 
for the two remaining governing equations at the boundary. This approach is used 
at a boundary surface (e.g. £ = £ max in Figs, la, lb) The second procedure is to 
specify two combinations of the dependent variables and use one sided differencing 
for a third governing equation at the boundary, this approach is used at an axis 
(e.g. £ = 0 in Fig. la). In both cases, periodicity is assumed in the 77 direction. 

The surface grid can be arbitrarily specified by the user. For relatively simple 
body shapes, the user can generate an orthogonal grid along the body surface. This 
process assures that the resulting grid will be nearly orthogonal in all directions. 
For complicated configurations, however, the surface grid is generated along cross 
sections of constant axial distance. Depending on the body shape, the £ and 77 
lines can be far from being orthogonal. However, this does not introduce any 
inconsistency because the orthogonality between the £ and 77 is not enforced in 
the grid generation procedure. Experience has shown that the overall quality of 
the mesh is dependent upon the quality of the surface grid. Therefore, the user 
should pay extra attention in defining the surface grid to ensure a well behaved grid 
everywhere. 

Since the grid is obtained by marching outward from the body, an exact location 
for the outer boundary can not be prescribed. Nevertheless, an approximate location 
for the outer boundary can be specified through proper selection of the cell volumes 
used in Eq. (lc). The cell volumes are defined by connecting straight lines between 
the surface boundary and the desired outer boundary shape. A prescribed clustering 
is then used to define the points along each f line. These cell volumes are then used 
in Eq.(lc) for the marching process. The resulting outer boundary will be close to 
the desired boundary. If a more exact outer boundary shape is needed, the computed 
grid points can be reclustered in the f direction to allow that. This is accomplished 
by calculating the intersection of the computed f lines with the desired boundary. 
The grid points are then obtained by redistributing points between the body and 


4 


the calculated intersection points along each line using the prescribed clustering. 
Although, the grid remains body normal, a significant amount of mesh skewness 
can be introduced especially if the desired outer boundary is relatively close to the 
body surface. This can be remedied by adding an extra global marching pass to 
improve the quality of the grid. This is done by computing new cell volumes at the 
end of the first global marching step and using that as the specified cell volume for 
the second global step. 

3. FLOW SIMULATION PROCEDURE 

Although in the current study we are dealing with predominantly supersonic 
flows, some regions of the flow field may contain pockets of embedded subsonic flow 
and/or axial separation. In such regions, the flow equations are elliptic, precluding 
the use of the more efficient space-marching techniques such as the parabolized 
Navier-Stokes (PNS) codes (e.g. Ref. [9]). Instead, a time-dependent procedure 
(F3D)is used to find the steady-state solution of the three-dimensional conservation- 
law form of the Reynolds-averaged Navier-Stokes equations. 

3.1 Numerical Procedure 

The unsteady Navier-Stokes equations with the thin-layer approximation in a 
strong conservative form, can be written as: 

Qt + + F n + G j = Re 1 (4) 

where Q is the dependent variables vector, E, F and G are the inviscid-flux vectors, 
and the S is contribution from the viscous terms. These vectors are given in Ref. [2]. 

A first or second order upwind differencing is used in f ( streamwise) direction 
while a second order central differencing is used in r) (meridional) and f (normal) 
directions. The implicit non-iterative two-factored scheme developed in Ref. [2] is 
used to solve Eq.(4). In this scheme, the flow variables at time step n + 1 are 
obtained from : 

£ + £~{Q n+1 - Q n ) = -At{6*{E+) n + 6({E~) n + 6'F n + 6'G n - Re-'S^S") (5) 

where, 6 C is the central operator, is the backward operator ,6* is the forward 
operator and 6 is the midpoint central operator. In Eq.(5), £ + is a lower block 
triangular matrix and £~ is an upper block triangular matrix and they are given 
by: 

£ + = (/+ hAt6^A + + 6^C - Re~ 1 6^{J~ 1 MJ)) 

£' = (/ + hAt6{A~ + 6‘B) 
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where the 5x5 matrices A ± ,B,C and M are the flux Jacobians of the fluxes 
E ± ,F,G and S respectively. 

Equation (5) is solved by a sequence of two one-dimensional like inversions. In 
practice, smoothing terms are added in the r) and f directions where central spatial 
differencing is used. The smoothing consists of second order implicit terms and 
combination of second and fourh order explicit terms. Also, free stream is subtracted 
as a base solution to reduce the numerical differencing error. 

In the above formulation, upwind differencing is used in the streamwise direction. 
In hypersonic flows, the bow shock is typically more aligned with f lines than with 
the £ lines especially near the nose section where the bow shock is the strongest. 
Therefore, it would be more advantageous to use upwind differencing in the f direc- 
tion and central differencing in the £ and 77 directions. Without making significant 
changes to the code, upwind differencing can be utilized in the outward direction 
by simply interchanging £ with f and neglecting the implicit portion of the viscous 
terms in Eq.(5) resulting in: 


£ + £~(Q n+ 1 - Q n ) = - A t{6l'{G + ) n + 6f{G~) n + 6‘F n +6lE n - Re- l S^S n ) (6) 

where, 

£+ = (/ + hAt6*C + + S^A) 
i~ = (I + hM6(C~ + 6‘B) 

Here C ± is the flux jacobian of G ± . It must be noted that, while upwind differencing 
is used for the inviscid flux in the f direction, the viscous flux on the right hand 
side of Eq.(6) is still centrally differenced using the midpoint operator. 

The above alternative formulation produces sharper bow shocks as will be shown 
later. Moreover, by interchanging the algorithm this way, the streamwise axis sin- 
gularity(at which the grid transformation Jacobian is infinite)is tied to central dif- 
ferencing rather than the most sensitive upwind differencing. The obvious drawback 
is that the viscous terms are treated explicitly. This proved to reduce the maximum 
allowable time step increment resulting in slower convergence rates. The possibility 
of modifyng the code in order to retain the implicit viscous terms in Eq.(6) will be 
addressed in a forthcoming report. 

The normal flux split formulation is implemented into the code as an option. An 
input control parameter activates the calls to the appropriate boundary conditions 
subroutine. In order to compare the two formulations, the inviscid flow around the 
shuttle nose at Mach number of =7.9 and 25 degrees angle of attack was selected 
as a test case. The grid in the pitch plane of symmetry is shown in Fig. 2. Figures 
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3 and 4 show the computed steady state mach number contours in the pitch plane 
of symmetry obtained using streamwise flux split (Eq.(5)) and normal flux split 
(Eq.(6)), respectively. It is clearly seen that the normal flux split yields a much 
sharper bow shock definition near the nose and in the windward region where the 
shock is roughly aligned with £=constant surfaces. However, both the streamwise 
and normal flux split result in shock smearing in the downstream portion of the 
leeside. This is because the bow shock there is no longer aligned with either the g 
or £ surfaces as could be seen from the grid in Fig. 2. It is worth noting that, with 
the normal flux split, a simple grid adaptation where the grid points are allowed 
to move on the f lines can solve this problem. In general, a formulation in which 
flux splitting is utilized in both the streamwise and normal directions may be more 
suitable for general hypersonic flow problems. 

3.2 Boundary Conditions 

At f = 0 (body surface), the viscous no-slip condition is imposed. The normal 
momentum equation is used to compute the surface pressure. In addition, either an 
adiabatic wall condition or a prescribed wall temperature are specified. At f = f maz 
(outer boundary), the bow shock can either be fitted or captured. In the bow-shock- 
capturing calculations, the boundary values are fixed at its free stream values. With 
the bow shock fitted, its location as well as the flow variables behind it are deter- 
mined from the shock fitting procedure of Ref. [10]. The shock is allowed to move 
along the curved f lines. The original clustering is maintained by redistributing 
points along f lines after each time step to adjust for the bow shock movement. 

In the present study only cases with bilateral symmetry are considered, and the 
boundary conditions in the meridional direction rj enforces symmetry conditions 
at the rj — 0 (leeward) and tj = r} max (windward) planes. The symmetry is im- 
posed through the use of reflection points on the opposite sides of the leeward and 
windward planes. 

At f = 0, either axis (”C-0” mesh) or inflow (”H-0” mesh) boundary conditions 
are imposed. At the axis the boundary conditions are obtained by extrapolating 
and then averaging the flow variables from the surrounding points. When upwind 
differencing is used in the streamwise direction, the flow variables at the first set 
of grid points away from the axis is obtained by interpolation to avoid integrating 
the equations of motion at these points. This interpolation step is not needed when 
central differencing is used in the streamwise direction. At the inflow boundary, all 
the flow variables are specified. At £ = f mai (outflow boundary) an extrapolation 
formula consistent with supersonic outflow is used. 

All the boundary conditions are implemented explicitly in the code. While the use 
of explicit boundary conditions with an implicit scheme does cause some degradation 
of the convergence rates, the use of explicit boundary conditions allows the user to 


7 



treat different flow problems by merely changing the boundary conditions routines 
without changing the rest of the code. 


4. RESULTS 


The current grid generation and flow solution procedures have been used to sim- 
ulate the hypersonic flow field around the shuttle orbiter and a generic wing-body 
configuration. 

4.1 Shuttle Orbiter 

Figure 5a shows an overall view of the mesh generated for the orbiter while Fig. 
5b shows the details of the grid near the body surface. The three dimensional mesh 
is displayed by showing the surface grid (f = 0), the grid in the plane of symmetry 
(n = 0 and y = r] max ) and the grid at the exit plane ( £ = £ mai ). It is seen that 
the grid is well behaved and is body normal. In generating the mesh, an estimate 
for the bow shock location was made and the appropriate values for the cell volume 
were chosen to result in an outer boundary location suitable for shock capturing 
computations. The boundary procedure corresponding to Fig. la (”C-0” mesh)is 
used where the two relations y = y(x) and z = z(x ) are employed at the axis. At the 
outflow boundary the relation x — constant is employed. The value of the constant 
is determined from the surface grid. The grid is compromised of 54 points along 
the body, 61 points in the cross section and 45 points in the outward direction. 

The flow conditions are = 7.9 and a — 25 degrees, and a constant wall 
temperature of 540 degrees. The Reynolds number (based on the orbiter total 
length) is 1.46x10°. The flow was assumed to be laminar. These flow conditions 
correspond to the wind tunnel experimental data of Ref. [11]. The calculations were 
performed using the streamwise flux split option. The result of the computations in 
the region where the wing shock interacts with the bow shock are shown in Figs. 6 
and 7. Figure 6 shows the pressure contours at a typical cross section. The pressure 
and mach number contours in a y = constant plane passing by the wing tip are 
shown in Fig. 7. The general features of the flow field near the surface is indicated 
in Fig. 8 which depicts the simulated oilflow pattern. Lines of cross flow separation 
are clearly visible on the upper wing and on the top part of the orbiter. Also, a line 
of axial flow separation is observed on the top surface. Fig. 9 show the variation of 
the surface pressure along the windward plane of symmetry and in a cross plane at 
x jl — .796. The experimental data at xjl — .8 is also shown, and it seen that there 
is a good agreement between the computed results and the data. The steady-state 
computations required about 7 MW of core memory and took about 8 hours of 
CRAY 2 time. 
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4.2 Generic Wing-Body Geometry 

The flow about a generic wing body configuration was also computed. The flow 
conditions for this case are M 0 0 = 25, a = 5 degrees, and a constant wall temper- 
ature of 3000 degrees. A specific heat ratio of 1.2 was chosen to simulate the real 
gas effects at high freestream Mach number. 

In order to facilitate the solution at this high Mach number, a segmentation pro- 
cedure(Ref.jl2]) in which the flow field is divided into two regions is used. The 
boundary between the two adjacent regions is carefully located beyond the rounded 
nose where the outer inviscid flow is supersonic and the boundary layer is attached. 
This ensures the absence of any upstream influence between the two regions. This di- 
vision enables the use of shock fitting in the nose region where the shock is strongest, 
while shock capturing is used in the downstream section. The use of shock captur- 
ing in the second segment avoids the tedious process associated with guessing the 
initial bow shock shape for the shock fitting procedure. 

The grid in the second segment was first created with the outer boundary located 
far away from the body. Then, the local bow shock slope at the end of first seg- 
ment was used as guide for determining a more appropriate location for the outer 
boundary. The grid points were then redistributed in the £ direction to reflect this 
new location as explained earlier. The boundary procedure corresponding to Fig. 
lb (”H-0” mesh) is used where the relation x — constant are employed at the inflow 
and the outflow boundaries. The surface and resulting computational grid (at the 
symmetry plane and exit plane) are shown in Fig. 10. The grid consisted of 45 
points in the streamwise direction, 92 points in the meridional direction and only 
26 points in the outward direction. It is seen that the grid is smooth and body 
normal everywhere except at £ = £ mQI where one of the orthogonality equations 
was replaced by the relationship x = constant to ensure stable marching. 

The computation was first performed with the original flux split formulation 
until steady state was reached. Then, the normal flux split option was activated to 
produce a better defintion for the bow shock. The implicit viscous terms on the left 
hand side of Eq.(5) were neglected as mentioned before. This forced a reduction 
in the time step by an order of magnitude thus slowing down the convergence 
considerably. Figure 11 show the Mach number in a cross plane near the end of the 
vehicle. The pressure contours on the body surface are shown in Fig. 12. It is seen 
that there is a pressure rise near the nose and near the wing tip. 

CONCLUDING REMARKS 

The hyperbolic grid generation scheme allows the generation of a three dimen- 
sional body normal mesh suitable for viscous flow computations. The limited outer- 
boundary control is suitable for hypersonic applications. A method for solving the 
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viscous hypersonic flow field around realistic configurations on the CRAY2 computer 
was presented. Flux splitting was used in the streamwise and normal directions. 
The normal flux splitting results in a sharper bow shock and avoids the axis singu- 
larity problem but requires more iterations because of the explicit treatment of the 
viscous terms. 
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(a) C-0 MESH (b) H-0 MESH 

Figure 1. Mesh topologies for hypersonic flow computations. 
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Figure 3. Mach number contou 
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Figure 9a. Comparison between the computed pressure and the experimental 
data on the windward side. 
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Figure 9b. Comparison between the computed pressure and the experimental 
data at x jl — 0.8. 
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Figure 11. Mach number contours in a £ constant plane near the end of the 
vehicle. 
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